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Abstract 

The definition and computation of the topological susceptibility in non-abelian gauge theo- 
ries is complicated by the presence of non-integrable short-distance singularities. Recently, 
alternative representations of the susceptibility were discovered, which are singularity-free 
and do not require renormalization. Such an expression is here studied quantitatively, using 
the lattice formulation of the SU(3) gauge theory and numerical simulations. The results 
confirm the expected scaling of the susceptibility with respect to the lattice spacing and 
they also agree, within errors, with computations of the susceptibility based on the use of 
a chiral lattice Dirac operator. 



1. Introduction 

In QCD and other non-abelian gauge theories, the discussion of the effects of the 
topological properties of the classical field space tends to be conceptually non-trivial, 
because the gauge field integrated over in the functional integral is, with probability 
1, nowhere continuous. The topological susceptibility, for example, is only formally 
given by the two-point function of the topological density at zero momentum, unless 
a prescription is supplied of how exactly the non-integrable short-distance singularity 
of the two-point function is to be treated. 

In lattice gauge theory, the problem was reexamined some time ago [1-3] starting 
from a formulation of lattice QCD which preserves chiral symmetry. An important 
result of this work was that the topological susceptibility can be written as a ratio of 
expectation values of other observables which remains well-defined in the continuum 
limit. A particular choice of regularization is then not required, i.e. the new formula 
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provides a universal definition of the susceptibility. Moreover, this definition is such 
that the anomalous chiral Ward identities are fully respected. 

The aim in the present paper is to complement these theoretical developments by 
demonstrating the suitability of the universal definition for the computation of the 
topological susceptibility in lattice gauge theory. In this study, the pure SU(3) gauge 
theory is considered and a recently proposed version [4] of the universal formula is 
used (see sect. 2). As far as the feasibility of the calculation is concerned, the results 
are however expected to be directly relevant for QCD too. 



2. Singularity-free expressions for the topological susceptibility 

The formula for the susceptibility obtained in [4] is not very complicated, but some 
preparation is required to be able to write it down. Prom the beginning, the theory 
is considered on a finite hypercubic lattice with spacing a, volume V and periodic 
boundary conditions. While some particular choices have to be made along the way, 
these details are expected to be irrelevant in the continuum limit in view of the fact 
that the expression is renormalized and free of short-distance singularities. 

2. 1 Spectral-projector formula 

The construction starts by adding a multiplet of valence quarks with bare mass mo 
to the theory. On the lattice, the added fields are taken to be of the Wilson type [5] 
and the associated massive Dirac operator is assumed to include the Pauli term 
required for 0(a) improvement [6,7] (the relevant improvement and renormalization 
constants are collected in appendix A). 

The hermitian operator Dm has a complete set of orthonormal eigenmodes 
with non-negative eigenvalues a. On average there are only few eigenvalues below 
some threshold cith proportional to the square of the valence-quark mass (see fig. 1). 
Above the threshold, the spectrum has an approximately constant density with a 
slight downward trend in the range considered in the figure. Siich a trend is absent 
in two-flavour QCD [4], but is qualitatively in line with the behaviour of the spectral 
density at next-to- leading order of quenched chiral perturbation theory [8]. 

The topological susceptibility is now given by [4] 

^ (TV{Fm}) (Tr{75FM}Tr{75PM}) , 
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Fig. 1. Average number of eigenvalues of D„i Dm below on a 64 x 32^ lattice 
with spacing a ~ 0.07 fm, plotted as a function of the renormalized value Mr of M. 
The renormalized valence-quark mass mn is about 25 MeV in this example. 

where Pm denotes the orthogonal projector to the subspace spanned by the eigen- 
modes of D^^ Dm with eigenvalues a < M^. It is taken for granted in this formula 
that is above the effective threshold ath of the spectrum and that the renormal- 
ized valence-quark mass mp as well as the renormalized value Mr of M are held 
fixed when the lattice spacing is taken to zero (cf. appendix A). 

2.2 Alternative expressions 

Equation (2.1) derives from a study of the renormalization and symmetry properties 

of the n-point correlation functions of the scalar and pseudo-scalar densities of the 
valence quarks. There exist various representations of the topological susceptibility 
of a similar kind, all having the same continuum limit. In particular, 
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(TV{Ml,}) (TV{75M2^}TV{75 
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(TV{75Ml^75Ml^}) 



'^^^ +0(a2) 



(2.2) 



for any function Mjv^ of Dj^Dj^ which is equal to unity in the vicinity of the spectral 
threshold and rapidly decaying above . The shape of the function can otherwise 
be chosen arbitrarily and only affects the size of the O(a^) corrections. 

In the numerical work reported later, M^^f is set to the rational approximation of 
the projector Pm previously used in [4] for the computation of the mode number in 



3 



Table 1. Ladiro iiaramctcrs. stalisiics and bare mass parameters 



Lattice 


/3 


a [fm] 


-^cnfg 


K 


aM 


48 X 24^ 


5.96 


0.0999(4) 


100 


0.134433 


0.03188 


64 X 32^ 


6.17 


0.0710(3) 


100 


0.135540 


0.02232 


96 X 48^ 


6.42 


0.0498(3) 


100 


0.135561 


0.01566 



two-flavour QCD. While this function is not exactly equal to unity in the vicinity 
of the spectral threshold, the effect of the deviation on the calculated values of the 
topological susceptibility is totally negligible with respect to the statistical errors. 
For the reader's convenience, the function is given explicitly in appendix B. 



3. Numerical studies 

The expression on the right of eq. (2.2) is a ratio of well-defined expectation val- 
ues that can in principle be computed through numerical simulations. In practice, 
the traces Tr{. . .} can normally not be evaluated exactly, but as explained in sub- 
sect. 3.3, they can be estimated stochastically with a moderate computational effort 
and without compromising the correctness of the final results. 

3. 1 Simulation parameters 

The studies reported in this paper are based on simulations of the lattice theory at 
three values of the inverse bare gauge coupling /3 = Q/qq (see table 1). A well-known 
deficit of all currently available simulation algorithms for non-abelian gauge theories 
(including the link-update algorithms used here) is the fact that the integrated auto- 
correlation times of quantities related to the topological charge are rapidly growing 
when the lattice spacing decreases [9,10]. In order to guarantee the statistical inde- 
pendence of the A'^cnfg gauge-field configurations used for the "measurement" of the 
topological susceptibility, the distance in simulation time of subsequent configura- 
tions was required to be at least 10 times larger than the relevant autocorrelation 
times. 

Physical units are defined through the Sommer reference scale rg = 0.5 fm [11]. In 
the range of the gauge coupling covered here, the conversion factor ro/a from lattice 
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to physical units was accurately determined by Guagnelli et al. [12]. The spacings 
of the three lattices thus decrease from roug hly 0.1 to 0.05 fm by factors of l/\/2, 
while the lattice sizes in physical units are approximately constant. 



3.2 Spectral projector parameters 

As already mentioned, the operator Wm is taken to be a rational approximation to 
the projector Fm- It thus depends on the valence-quark mass, the mass M and the 
parameters n and e that determine the accuracy of the approximation (cf. appendix 
B). A reasonable choice of the latter, previously made ref. [4], is n = 32 and e = 0.01. 
In the range of eigenvalues of D^J Dm below 0.85 x M^, the approximation error is 
then smaller than 2.2 x 10~^, which is by far small enough to guarantee the absence 
of significant systematic effects in eq. (2.2). Moreover, the contribution of the high 
modes is safely suppressed. 

The valence-quark mass and the mass parameter M were adjusted such that their 
renormalized values in the MS scheme at normalization scale fi = 2 GeV are about 
25 and 100 MeV, respectively. Using the information collected in appendix A, the 
corresponding values of the bare mass parameters, k = {8 + 2amo)~^ and aM, can 
be worked out and are listed in table 1. On the lattices considered, there are then 
57 — 70 eigenmodes of D^J Dm with eigenvalues below and an average density 
of roughly 1 such mode per fm^. 

As already emphasized, the calculated values of the topological susceptibility are 
not expected to strongly depend on all these details and should in any case always 
extrapolate to the same value in the continuum limit. The lattice effects are unlikely 
to be small, however, if aM is not much smaller than 1 or if the expectation values 
on the right of eq. (2.2) would be dominated by the modes up to and slightly above 
the spectral threshold, where the effects arc kincmatically enhanced. Both of these 
unfavourable situations are avoided by the above choice of the mass parameters. 

3.3 Random-field representation 

In lattice QCD, random field representations were introduced many years ago [13] 
and are now widely used. The application of the method in the present context 
requires a set r/i, . . . , t/at of AT pseudo-fermion fields to be added to the theory with 
action 



where the bracket {rj, Q denotes the obvious scalar product of such fields. For every 



N 




(3.1) 



k=l 
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gauge field configuration, these fields are generated randomly so that one obtains a 
representative ensemble of fields for the complete theory. In the rest of this section, 
expectation values are always taken with respect to both the gauge field and pseudo- 
fermion fields. 

The stochastic observables 

1 ^ 

fe=l 
1 ^ 

H = — ^ (MMTs^MT/fe, Mm75Mm%) , (3.3) 
fe=i 

1 ^ 

C = - XI (Km%, 75Km%) , (3.4) 
fc=i 

may now be introduced and a moment of thought then shows that the expectation 
values on the right of eq. (2.2) are given by 

{Tr{RU) = (A), (3.5) 

(TV{75ML}T¥{75M2^}) = - i^, (3.6) 

(Tr{75ML75Mi,}) = {B). (3.7) 



The topological sTisccptibility can thus be computed by calculating the expectation 
values of A, B and C^. For a given gauge- field configuration, the evaluation of these 
observables requires the fields ^uVk, ^j^Vk and ^M^b^MVk to be computed, i.e. the 
total numerical effort per configuration is roughly equivalent the one required for 
3N applications of the operator Rm to a given pseudo-fermion field. 

From this point of view, small values of A'^ are favoured, but a good choice of N 
must also take into account the fact that the variance of the stochastic observables 
decreases with N. Some experimenting then shows that setting = 6 is a reasonable 
compromise at the specified values of the mass parameters. Since Rm is a rational 
function of D^D^ of degree [2n + 1, 2n + 1], the measurement of the stochastic 
observables requires the (twisted-mass) Dirac equation to be solved for altogether 
2340 source fields. The computational load thus tends to be heavy, but the problem 
is well suited for the application of highly efficient solver techniques such as local 
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Table 2. Mode number and topological susceplibilKy 



/3 




x\" [MeV] 


(xy')wF [McV] 




5.96 


1.053(18) 


197.3(7.7) 


187.7(6.0) 


1.051(25) 


6.17 


1.075(22) 


203.8(9.1) 


192.8(7.0) 


1.057(22) 


6.42 


1.060(24) 


186.6(9.9) 


181.0(7.3) 


1.031(26) 



deflation [14] (see ref. [15] for a recent review of the subject). In particular, when 
these are used, the effort scales linearly with the lattice size and is nearly independent 
of the values of the mass parameters. 

5.^ Simulation results 

The simulation data discussed in the following paragraphs are summarized in table 2. 
In all cases, the statistical errors were estimated using the jackknife method and were 
combined in quadrature with the quoted scale errors (where appropriate). 

(a) Mode number. The average number u of eigenmodes of Dr,^Dm with eigenvalues 
below is an extensive quantity and is therefore normalized by the lattice volume 
in table 2. At the specified bare masses, the renormalized masses and Mr are 
practically equal to 25 and 100 MeV, respectively, on all three lattices considered. 
In view of the renormalization properties of the mode number [4], the calculated 
values of u/V are thus expected to be the same up to O(a^) effects. 

Within errors, the figures listed in the second column of table 2 in fact coincide 
with one another. Note that the quoted errors do not take into account the fact that 
the mass renormalization factors and thus the renormalized values of the masses are 
only known up to an error of about 2% (see appendix A). Once this error is included 
in the analysis, one can still conclude, however, that the simulation results confirm 
the expected scaling of the mode number to the continuum limit at a level of precision 
of 3% or so. 

(b) Topological susceptibility. The values of the susceptibility calculated along the 
lines of the present paper are listed in the third column of table 2. Again one observes 
no statistically significant dependence on the lattice spacing. Finite-volume effects 
are, incidentally, known to be negligible with respect to the statistical errors on all 
three lattices [16]. 

Fits of the data by a constant or a linear function in yield consistent results 
in the continuum limit. Since the slope in turns out to vanish within errors, the 
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(more accurate) number 

= 196.5(5.1) MeV (3.8) 

obtained by fitting with a constant is quoted here. This result happens to be prac- 
ticahy on top of the value 194.5(2.4) MeV obtained by Del Debbio, Giusti and Pica 
[16] using a chiral lattice Dirac operator and the index theorem [17] f. 
(c) Charge sectors & the Wilson flow. An understanding of how exactly the topolog- 
ical charge sectors emerge in the continuum limit has recently been achieved using 
the Wilson flow [18]. The definition of the topological susceptibility suggested by the 
sector division is geometrically appealing and computationally far less demanding 
than the spectral-projector formula (2.2). Presumably the two definitions agree in 
the continuum limit, but there is currently no solid theoretical argument that would 
show this to be the case. 

The values of the susceptibility computed using the Wilson flow are listed in the 
fourth column of table 2 (see ref. [18] for the details of the calculation). While they 
appear to be systematically lower than the ones obtained using the spectral-projector 
formula, the differences are statistically insignificant on each lattice. Moreover, there 
could be lattice effects of size up to the level of the statistical errors. 

Since the same ensemble of representative gauge-field configurations was used in 
the two cases, the quoted errors are correlated to some extent (not completely so 
in view of the use of random fields). The ratio listed in the last column of table 2 
is therefore obtained with slightly better precision than the susceptibilities. Fits of 
the ratio by a constant and linear function in arc both possible, the values in the 
continuum limit being 1.048(14) and 1.036(31), respectively. The spectral-projector 
and the Wilson-flow definition of the susceptibility thus coincide to a precision of a 
few percent. While there is some tension in the data, there is no clear evidence for 
the definitions to be different at this level of accuracy. 



4. Conclusions 

The numerical studies reported in this paper confirm the universality of the spectral- 
projector formula (2.2) for the topological susceptibility. In particular, no statisti- 
cally significant lattice-spacing effects were observed and the calculated values agree 

f In ref. [16], a different convention for the conversion from lattice to physical units was used and 
the value for the susceptibility quoted there is therefore slightly different from the one printed here. 
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with the result obtained by Del Debbio, Giusti and Pica [16], where a chiral lattice 
Dirac operator was used. 

The numerical effort required for the stochastic evaluation of the spectral-projector 
formula increases proportionally to the number V ja^ of lattice points, but is a flat 
function of all other parameters. On large lattices, computations of the susceptibility 
along these lines thus tend to be more feasible than those based on a chiral lattice 
Dirac operator (which scale roughly like V^). Even less computer time is required if 
the susceptibility is defined via the Wilson flow, but a formal proof for this definition 
to be in the same universality class as the spectral-projector formula is still missing. 

With respect to the pure gauge theory, the application of the spectral-projector 
formula in QCD is not expected to run into additional difficulties. Accurate cal- 
culations of the topological susceptibility however require representative ensembles 
of, say, a few hundred statistically independent gauge-field configurations to be gen- 
erated. This part of the calculation usually consumes most of the computer time 
and may rapidly become prohibitively expensive at small lattice spacings [10]. At 
present, computations of the susceptibility on lattices similar to the ones considered 
here are therefore not easily extended to QCD with light sea quarks. 

We wish to thank Leonardo Giusti for helpful discussions on various issues related 
to this work. All numerical calculations were performed on a dedicated PC cluster 
at CERN. We are grateful to the GERN management for providing the required 
funds and to the CERN IT Department for technical support. F. P. acknowledges 
financial support by an EIF Marie Curie fellowship of the European Community's 
Seventh Framework Programme under contract number PIEF-GA-2009-235904. 



Appendix A. 0(a) improvement and renormalization 

A.l Dirac operator and renormalization constants 

The lattice theory considered in this paper is set up as usual, using the Wilson gauge 
action and the standard 0(o)-improved Wilson-Dirac operator. The notation and 
normalizations are as in ref. [7]. In particular, and ca denote the coefficients of 
the Pauli term in the Dirac operator and the 0(a) term required for the improvement 
of the axial current. Here these coefficients were set to the values given by the non- 
perturbatively determined interpolation formula quoted in ref. [19] (see table 3). 

The values of the renormalization constant Za of the axial current listed in table 3 
were obtained by evaluating the interpolation formula given in ref. [20]. In the case 
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Table 3. Improvement coefHcients and renormalization constants 







CA 


Za 


Zp 


5.96 


1.81663 


-0.11432 


0.789(8) 


0.629(14) 


6.17 


1.63125 


-0.04015 


0.807(8) 


0.622(14) 


6.42 


1.51877 


-0.02446 


0.824(8) 


0.618(13) 



of the renormalization constant Zp of the pseudo-scalar quark density, the quoted 
values are the ones required to pass from the lattice normalization of the density to 

the one in the MS scheme of dimensional regularization at normalization scale fi = 2 
GeV. The constant was calculated in two steps, first passing from the lattice to the 
renormalization- group- invariant normalization [21] and then from there to the MS 
scheme [22]. 

A. 2 Quark masses 

The renormalized quark mass in the MS scheme is given by 

Zp{l + bpam^)"'^^^''>' ^^-^^ 

where m is the bare current-quark mass, rriq = mo — rUc the subtracted bare mass 
and rric the critical bare mass. Here and below, bx (where X = A,P,...) denotes 
an improvement coefficient required to cancel lattice effects proportional to aniq. 

At fixed gauge coupling, the current quark mass is related to the subtracted bare 
mass through 

m = Zmq{l + (6™ -bA + bp)am^} + O(a^). (A.2) 

Using the Schrodinger functional, the coefficients Z and bm — bA+bp were determined 
non-perturbatively by Guagnelli et al. [23] (see table 4). Also shown in tabic 4 is 
the current quark mass at one value of the hopping parameter k = {8 + 2amo)~^. 
Together with eq. (A.2), these data allow the current quark mass to be estimated 
at larger values of k, where a direct computation on the lattices considered in this 
paper tends to be compromised by the presence of accidental near-zero modes of the 
Dirac operator. 
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Table 4. Quark mass reference point and extrapolation coefficients 





K 


am 


Z 


hm-hA + hp 


5.96 


0.13360 


0.033216(60) 


1.0402(4) 


-1.017(13) 


6.17 


0.13521 


0.016636(39) 


1.0935(4) 


-0.739(13) 


6.42 


0.13545 


0.008046(20) 


1.1041(4) 


-0.687(13) 



A. 3 Renormalization of the mode number 

The renormalization and improvement properties of the mode number 

KM,mq) = TV{Pm} (A.3) 

were discussed in detail in ref. [4] . In particular, it was shown there that 

z/r(Mr, mR) = iy{M, m^) (A.4) 

is a renormahzed and 0(a)-improved quantity. In this equation, the bare parameters 
M, rriq are to be expressed through the renormahzed masses 

Mr = Zp\l + bi,amq)M (A.5) 

and rriR. Note that Mr does not renormahze in the same way as the quark mass. 
Currently the coefficient 

6^ = -i- 0.111(4) X 32 + 0(^4), (A.6) 
is only known to one-loop order of perturbation theory [24,4]. 



Appendix B. Definition of Rm 

The operator M.m is of the form [4] 

RM = hiX), X = l- - (B.l) 
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where M* 2± M and h{x) is a polynomial approximation to the step function 0{—x) 
in the range — 1 < x < 1. More precisely, the polynomial is given by 

h{x) = \{l-xP{x^)] , (B.2) 

P{y) being the polynomial of degree n which minimizes the deviation 

5= mEQC |l-v^P(y)| (B.3) 

e<y<l 

for some specified (small) value of e. This choice ensures that h{x) provides a uniform 
approximation to the step function in the range > y^, with maximal absolute 
deviation equal to \5. Moreover, inspection shows that h{x) decreases monotonically 
in the transition region —^/e<x< \/e. 

For a given degree n and transition range e, the coefficients of the minmax poly- 
nomial P{y) can be computed mimcrically using standard techniques. An efficient 
procedure was described in ref. [25], for example. The mass oc M is then deter- 
mined through 

w. - (^) + - ' + °(^)- 

As explained in appendix B of ref. [4], this convention is intended to minimize 
the deviation |Tr{PM — I^mII- the present context, other choices of M* would 
however do just as well, since eq. (2.2) is expected to hold for any M. 

Small approximation errors 5 are achieved with moderately high degrees if e is 
not very small. For n = 32 and e = 0.01, for example, one obtains 

5 = 4.37 X IQ-^, M/M* = 0.96334. (B.5) 
The transition range |a;| < approximately corresponds to the range 

0.9 < Va/M* < 1.1 (B.6) 
of eigenvalues a of D^D^ in this case. 
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